Differential Equations

Introduction

Previously, we have seen many types of differential equations and have even solved some simple ones. However, what do we do if the equation is too hard for us to solve? What if we cannot express the answer in a convenient form? In these cases, mathematics software is an extremely useful tool.
In this notebook, we will use Mathematica to solve differential equations of the form (dy/dx)=f(x,y). We will learn how to express differential equations in a form that Mathematica can understand, as well as how to solve these equations symbolically. We will also learn how to represent the equation graphically, with direction fields. Finally, we will focus on initial value problems and how to solve them, both with direction fields and with the classic Euler's method.

Syntax

In order to take a derivative, we use the differential operator D. The first argument is the function of which to take the derivative, and the second argument is the variable with which the derivative is taken. This resembles Leibniz notation.

D[Sin[x],x]
D[x^2,x]
f[z_]:=Sqrt[z];D[f[z],z]

We can also use prime notation when we have a function of one variable:

Sin'[x]
f'[z]

Using prime notation can cause confusion when there are multiple variables present, so we are going to use D in this notebook. Differential equations are expressed the same way as algebraic equations, with two equals signs:

D[Sin[x],x]==0
D[x^2,x]==0
D[y[x],x]==f[x]

The Solve command does not work for these equations, unfortunately. Instead, we have to use the command DSolve:

Solve[D[y[x],x]==f[x],y[x]]
DSolve[D[y[x],x]==f[x],y[x],x]

With only the differential equation, this is as close as we can get to any form of solution. We don't have just one answer: we have many, many answers, and which answer is dependent on the value of the constant C. This statement of many answers and their dependence is called the General Solution of the differential equation. In order to narrow our answer down to one function of x, we need to solve a related problem called the initial value problem. An initial value problem looks like this:
	(dy/dx)=f(x,y), y(a)=b
Not only do we have the differential equation, but we also have that our answer must pass through a specific point (a,b). This point is called the initial condition or initial value. Initial conditions go with the original equation, in a list:

DSolve[{D[y[x],x]==f[x], y[1]==1},y[x],x]

As we can see, we now have only one solution, a function only of x.

Direction Fields

Sometimes the output of DSolve isn't useful. The equation is just too complicated.

f[x_,y_]:=(1-x^2)*y-x

DSolve[D[y[x],x]==f[x,y[x]],y[x],x]

In these situations, solving it by hand is usually out of the question as well. This is where numerical solutions come in. The first technique is plotting a special graph called a direction field or slope field. It indicates how the function changes, given the current state. Just as plotting the points helps us visualize the equation y=f(x), plotting the directions helps us visualize the equation (dy/dx)=f(x,y). We obtain this field using the command VectorPlot. VectorPlot is meant to visualize velocity fields. It takes three agruments:

The first argument is a list showing the velocities (or rates of change) in the x and y directions, respectively. For a differential equation of the form (dy/dx)=f(x,y), this will always be {1,f[x,y]} since (dx/dx)=1.

The second and third arguments are the range over which to plot. Just like with the Plot command, these are of the form {variable, lowerbound, upperbound}.

VectorPlot[{1,f[x,y]},{x,-3,3},{y,-3,3}]

By default, the lengths of the arrows are dependent of the relative magnitude of the slopes. This can make the plot very difficult to read since some arrows are very small.

In order to change the relative sizes of the arrows, we use the VectorScale option. VectorScale is passed a list: the first item in the list is the thickness of the arrow as a whole, the second is the thickness of the head only, and the third is how they scale with magnitude. Useful options are Tiny, Small, Medium, Large, and Automatic for the first two, and None and Automatic for the last.

VectorPlot[{1,f[x,y]},{x,-3,3},{y,-3,3},VectorScale-> {Automatic, Automatic, Automatic}]
VectorPlot[{1,f[x,y]},{x,-3,3},{y,-3,3},VectorScale-> {Tiny,Large,None}]

We can make more arrows with the VectorPoints option:

VectorPlot[{1,f[x,y]},{x,-3,3},{y,-3,3},VectorPoints-> 30,VectorScale-> {Tiny,Large,None}]

We can imagine a particular solution to an IVP as a curve which passes through the initial point and follows the direction of the arrows at all points along the curve. We can plot a few of these solutions by passing the StreamPoints option:

VectorPlot[{1,f[x,y]},{x,-3,3},{y,-3,3},StreamPoints-> 20]

Giving StreamPoints a list of points allows us to specify which points these curves must pass through:

VectorPlot[{1,f[x,y]},{x,-3,3},{y,-3,3},StreamPoints-> {{0,0}}]

Euler's Method

Let's return to the last sentence of the previous section: the solution to an initial value problem is a curve which passes through the initial point and follows the direction of the direction field at all the points along the curve. We can obtain a solution that obeys this rule by an iterative algorithm. Given our initial value problem, we can obtain a collection of points which represent the curve of our solution. Euler's method is one such way of obtaining these points, and is similar to the left-hand endpoints method of numerical integration. Euler's method is described by the following:

Given: the initial point ((Subscript[x, 1]),(Subscript[y, 1])), the differential equation (dy/dx)=f(x,y), and a parameter called the step size, h

Obtain the next point, ((Subscript[x, 2]),(Subscript[y, 2])) by:
	(Subscript[x, 2])=(Subscript[x, 1])+h
	(Subscript[y, 2])=(Subscript[y, 1])+f((Subscript[x, 1]),(Subscript[y, 1]))*h
This can be thought of as marching a small way along the tangent line to the solution curve at ((Subscript[x, 1]),(Subscript[y, 1])) until (Subscript[x, 2]) is reached, and using the value for the tangent line as the value for the solution curve.

Repeat step 2. to obtain ((Subscript[x, 3]),(Subscript[y, 3])), ((Subscript[x, 4]),(Subscript[y, 4])), etc.

For example, let's compute the solution to the equation (dy/dx)=(1-(x^2))y-x with the initial condition y(1)=1. We'll set h=(1/5), and our number of steps will be n=10: We are going to want x and y to be lists with n+1 (n steps plus the one initial condition) elements so that it will be easy to plot them later.

f[x_,y_]:=(1-x^2)*y-x
n=10
h=1/5
x=ConstantArray[0,n+1];
y=ConstantArray[0,n+1];
x[[1]]=1
y[[1]]=1

Now let's iterate over the steps. We will use the variable i to indicate the current step:

Do[
x[[i+1]]=x[[i]]+h;
y[[i+1]]=y[[i]]+f[x[[i]],y[[i]]]*h;,
{i,1,n}
]

We now have two lists: a list of the x values and a list of the y values. In order to plot them, we need a list of x,y pairs. We accomplish this with the Table command:

points=Table[{x[[i]],y[[i]]},{i,n+1}];

To plot them, we need to use the ListPlot command:

ListPlot[points]

The points we have can be connected using the Joined option:

ListPlot[points, Joined-> True]

A side effect of Euler's Method is that the accuracy is greatly dependent on the number of steps that you use. Let's try it again with 1000 steps:

f[x_,y_]:=(1-x^2)*y-x
n=1000;
h=1/500;
x=ConstantArray[0,n+1];
y=ConstantArray[0,n+1];
x[[1]]=1;
y[[1]]=1;
Do[
x[[i+1]]=x[[i]]+h;
y[[i+1]]=y[[i]]+f[x[[i]],y[[i]]]*h;,
{i,1,n}
]
points1000=Table[{x[[i]],y[[i]]},{i,n+1}];
ListPlot[points1000, Joined-> True]

Comparing the two answers:

ListPlot[{points,points1000},Joined-> True]

Mathematica can automate the process of solving Initial Value Problems with the NDSolve command:

Clear[x,y]
pointsND=y[x]/.NDSolve[{D[y[x],x]==f[x,y[x]],y[1]==1},y,{x,1,3}]
Plot[pointsND,{x,1,3}]

NDSolve does not use Euler's method by default. However, we can tell it to do so. We also have to specify the step size:

pointsNDeuler=y[x]/.NDSolve[{D[y[x],x]==f[x,y[x]],y[1]==1},y,{x,1,3},Method-> "ExplicitEuler",StartingStepSize-> 1/500]
Plot[pointsNDeuler,{x,1,3}]

Comparison with the answer we obtained via iteration:

Show[{ListPlot[{points1000},Joined-> True,PlotStyle-> Green],Plot[pointsNDeuler,{x,1,3},PlotStyle-> Red]}]

We can also compare with the direction field:

Show[{VectorPlot[{1,f[x,y]},{x,1,3},{y,-1,1},VectorStyle-> Red,VectorPoints-> 30,VectorScale-> {Tiny,Large,None}],Plot[pointsND,{x,1,3},PlotStyle-> Thickness[0.01]]}]

Exercises

Assignment 1: Using Direction Fields

Use Mathematica to plot the direction field of the differential equation (dy/dx)+y=6+cos(4x)-2sin(6x), and investigate the behavior of solutions to this equation. Your answer should include:

A direction field plot over the region -10≤x≤20,  -10≤y≤10 with all options at default.

A direction field plot over the region -10≤x≤20,  -10≤y≤10 with the VectorPoints option set to at least 100. In order for the arrows to be clearly visible, you may have to resize the plot window and/or set the scaling options.

A direction field plot over the region -10≤x≤20,  -10≤y≤10 with the StreamPoints option set so that at particular solutions passing through the following five points are visible: (0,-1), (0,0), (0,1), (0,2), and (0,3).

A conjecture based on the direction fields as to the behavior of the solution curves as x becomes large.

A general solution to the differential equation as produced by DSolve.

An explanation of why the general solution supports or contradicts the conjecture obtained from the direction fields.

Assignment 2: Using Euler's Method

Compute the solutions to the IVP (dy/dx)=(((e^x)-2y)/x), y(-2)=c for different values of c, and investigate how c affects the solution. Your answer should include:

Plots of the solutions obtained via Euler's method for c=0 and c=0.3 over the region -2≤x≤3 with a step size of at most 0.5.

A comparison between the solution curves for  c=0 and c=0.3.

An estimate of the critical point k at which the solution switches from being "c=0-like" to "c=0.3-like."

A plot of the solutions for k-0.01 and k+0.01 on the same graph. The two solutions should exhibit different behavior.

Assignment 3: The Importance of Step Size

Compute the solutions to the IVP (dy/dx)=(x^2)+(y^2), y(1)=1 using Euler's Method of different step sizes. Your answer should be a plot of all of the following on the same graph:

The solution via Euler's Method with h=0.1, n=10.

The solution via Euler's Method with h=0.02, n=50.

The default solution given by NDSolve over the same region.

Assignment 4: Comparing Direction Fields to Euler's Method

Solve the IVP (dy/dx)=y, y(0)=1 using both direction fields and Euler's Method, and compare the results. Your answer should include:

A direction field plot over the region -3≤x≤3,  -3≤y≤3, with the StreamPoints option set so that the solution to the IVP is also plotted.

A sketch of the solution curve, made by hand, by tracing over or writing on the direction field plot.

Plots of the solution curve over the same range, obtained via Euler's method with h=0.001 and h=0.0005.

An estimate of y(1) obtained for each of the three answers (both the Euler's method plots, and the graphically from the direction field plot).

A comparison of the three answers to the exact answer for y(1).

Assignment 5: Application to Population Growth

One of the ways of modeling the population of a species is with the logistic equation:
	(dP/dt)=kP(1-(P/K))
where P is the population, k is a constant called the growth factor, and K is a constant called the capacity. We can model a constant decrease in population, say from hunting, by adding an additional constant term, let's call it H:
	(dP/dt)=kP(1-(P/K))-H.
Produce direction fields for these two equations under the assumption that k=0.9, H=18, and K=1000, and compare the results. Your answer should include:

The two direction fields, plotted over a t range of zero to forty and a P range of -10 to 1100. You may need to resize the plot and/or adjust the scaling so that the arrows are clearly visible.

A comparison of the two direction fields.

An explanation of why the extra term H had this effect.